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ABSTRACT 



This thesis analyzes the nonlinear characteristics of motion stability of a 
submersible vehicle in combined sway, yaw, and roll motions. Previous results, at zero 
pitch angles, indicate that limit cycles are generated as a result of loss of stability. In this 
work, these results are extended to include nonzero pitch angles. This analysis can 
determine how changes in vehicle parameters and loading conditions will affect its 
operation and performance. Stability domains are generated for a variety of vehicle and 
environmental parameters. A nonlinear analysis is conducted in order to assess the 
stability characteristics of the resulting limit cycles. The results can lead to design 
guidelines for improving vehicle operational envelopes. 
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I. INTRODUCTION 



A. DESCRIPTION OF THE PROBLEM 

The study of dynamic responses in a submersible vehicle using a nonlinear 
analysis is important in determining operating envelopes for the vehicle. Previous studies 
have shown that in straight line motion, the coupling of the sway, yaw, and roll equations 
produce oscillating losses of stability in the system. [Ref. 5] Introducing a nonzero pitch 
angle to the vehicles motion will allow us to study the changes to the stability domain for 
a variety of environmental conditions. 

In our work, the linearized equations of motion are studied using an eigenvalue 
analysis to determine the systems stability through a variety of operational parameters. A 
nonlinear analysis is then conducted to assess the stability characteristics of the resulting 
limit cycles and their impact on the operating characteristics of the vehicle. 

B. OBJECTIVES AND OUTLINE 

In this thesis we expand on previous thesis work which examined the problem of 
stability in straight line motions of a submersible vehicle [Ref. 13]. The primary cause 
for this loss of stability is the coupling between the sway/yaw/roll equations of motion for 
a submersible vehicle. We know that the loss of stability creates stable limit cycles in 
straight line motion. This work analyzes the effects of introducing a nonzero pitch angle 
to the generic equations of motion in order to determine their effect on the creation of 
limit cycles. 
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The model used for this work is a variant of the Swimmer Delivery Model used in 
[Ref. 5] for which a generic set of hydrodynamic and geometric properties are available. 
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II. EQUATIONS OF MOTION 



A. COORDINATE SYSTEM 



A moving coordinate system was used for our analysis with the origin fixed at the 
vehicles center of buoyancy. The x-axis is fixed to the longitudinal plane of symmetry for 
the vehicle, the y-axis is positive to starboard, and the z-axis is positive downward. All 
symbols used in the development of the equations of motion are summarized in Table 1. 

B. GENERAL FORM OF THE EQUATIONS OF MOTION 

The equations of motion for a submersible vehicle in the horizontal plane are 
written as follows: 

Sway equation: 




'^pP + 1 "/ + Yp<, pq + y^rqr + + Y^vw + 

+ y,v + Y^Up + Y,Ur + Y,^vq + Y^^wp + Y^,wr + 
(W - B)cos dsiruf) - 




( 1 ) 
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Yaw equation: 



+ (^zz - ^yy W - ^xy (p‘ “ )" ^ yz (P'' + ^)+ K “ ^) + 

m[xc (v + f/r — wp)— (i/ - vr + w^)] = 

pq + A^^,^r + A^,t/v + N^vw + 

+ A^v^W + A^».p^P + + 

(xgW - j:g5)cose sin 0 + (y^W^ - >’fi5)sin 0+6^ - 

J^"°” /i(xXv + xr)^ + Co. b{x\w - xqf l^-^^xdx 
“■' ' ^cfW 



( 2 ) 



Roll equation: 

^xxP + (^zz - V + ^xy (p^ - ^)- hz )" ^xz (p^ + '") + 

mlydw-Uq + vp)- zdv + Ur - wp)] = 

K.p + K^r+Kp^pq + K^^qr + K^Uv + K^vw + 

K,ropU^ + K,v+K^Up + K^Ur + K,^vq + K^^wp + K„,wr + 

(y^W - y„5)cos0 cos0 - (zglY - z^R)cos 0 sin0 

(3) 

The rotational velocity equation around the x-axis: 

<i> = P (4) 

Ucf denotes the cross-flow velocity: 

U^^{x) = ^|{v + xrf+{w-xqf 
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Cn 


quadratic drag coefficient 


dr 


rudder deflection 


h(x) 


local height of hull 


(I XX yy zz) 


vehicle mass moments of inertia about body axes 


xyy^ yz*^ zx) 


cross products of inertia 


(KMN) 


moment components along the three axes 


m 


vehicle mass 


(p,q,r) 


rotational velocity components along the body axes 


(f>q.y) 


Euler angles 


u 


constant vehicle speed along the x -axis 


(u,v,w) 


translational velocities about (x.y.z) axes 


(x,y,z) 


distances along the three body axes 


(XXZ) 


force components along the body axes 


(xc<y G’^g) 


coordinates of the center of gravity 


(x B>y B’Zb) 


coordinates of the center of buoyancy 


^ nose 


fore coordinate of vehicle body 


^ tail 


aft coordinate of vehicle body 



Table 1: Nomenclature 



C. SIMPLIFICATIONS 

We must simplify the equations of motion in order to reflect the fact that we are 
analyzing motion about the y-axis. The simplifications that we employ are: 



• Acceleration, vv , in the z-direction is zero. 

• Acceleration in the longitudinal direction, u , is zero. 

• Rotational velocity, q, and acceleration, q , in the y- direction are zero 
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D. SIMPLIFIED EQUATIONS OF MOTION 



After applying the above simplifications, the equations of motion become: 



Sway equation: 

m\y + Ur-wp + Xcr-yc[p^ +r^)+Zap]= 

YpP + Y,r + Y,Uv+Y^vw+Y,U^5,+Y,v + Y^Up + Y^Ur + Y^^wp + Y^^wr + 
(W-5)cos6>sin^«>- 

[Cp h{xXv + xrf + Cp b{x)w^ ItT^ 

J-.au ^ - U^{X) 



Yaw equation: 

-IyzPr + KP + 

m[xc {v + Ur- wp)- ([/ - vr)]= 

N pp + N -r + N ^Uv+ N^vw+ Ng^ +N^v + N ^Up + N^Ur + wp + 

wr+{xcW - XgB)cos 0 sin 0 + (ygl^ - yB5)sin d + U ^ - 

f”” [Cd^^(.^Xv + ^0^ +Cob{x)w^^k^^xdx 

J-O.U > - U^f[x) 



Roll equation: 



^xxP + ^xyP’’-^yz^^ ~^xzj' + ^bc M~ Zq {v + Uk- Wp)] = 

Kpp + K/+K^Uv+K^vw+K^,„^U^+K,v + K^Up + K^Ur + 

K y^pWp + K ^^wr + [y^W - y g5)cos 6 cos (j>-{Z(;W- Zg R)cos G sin 0 
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Roll rate: 



(j) = p 



( 9 ) 
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i 




I 

i 
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III. LINEAR ANALYSIS 



A. LINEARIZATION 

The simplified equations of motion can be written in the matrix form: 

Ax = Bx + g(x) 



( 10 ) 



where the state vector, x , is defined as. 



and the state matrices are. 



X = 



A = 



m-y, 


mxc -Y- 


1 

1 


0 


mxc -N^ 




1 

1 


0 


-mZa-K, 






0 


0 


0 


0 


1 



and 



B = 



Y +Y w 
N ”h w 

V vw 

K.U + K^w 
0 



Y^U+Y^,-mU 
-mXfyU + N^U + N^^w 
mZcU + K^U + K^^w 
0 



YpU + Y^w+mw 0 

N^U + N^^w 0 

-mZrW+ K U + K^„w 0 



wp 



1 



0 
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The g(x) term remains as given [Ref. 13]. 

The nonlinear terms are linearized around a nominal point: 

^o=h'o^^O’Po^<t>oJ =0 



A Taylor series expansion is applied to the nonlinear terms about the nominal point, xq , 
and keeping only the linear components, the equations of motion, written in matrix form, 
become: 



Ax = B'x (^1) 

where, 

A' = A 





Y^U-mU 


YrfJ 


(W-B)cos0 


N^U 


-mx^U + N^U 


N^U 


(xgiy — XjR)cos0 


K^U 


mZaU + K^U 




• (-ZcW + ZgR)cos0 


0 


0 


1 


0 



To assess the dynamic stability of the vehicle, an Eigenvalue analysis is performed in the 
next section. 
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B. LOSS OF STABILITY 



An Eigenvalue analysis is used to determine the stability of the linearized system. 
Stability is dependent on the location of the Eigenvalues or the roots of the characteristic 
equation: 

det(5'->U' = 0) (12) 

in the polynomial form: 

AAV5A^+CA^+Z)A + £ = 0 (13) 

The coefficients equation (13) are given in [Ref. 13, pg. 11]. Using Routh’s criterion we 
can examine the stability of the system. The following conditions must be applied to the 
characteristic equation (13) to ensure all roots have negative real parts: 

BCD-AD^ -EB^ >0 (14) 

E>0 (15) 

If £■ is less than zero, one real root of (13) becomes positive and the system will become 
unstable in a divergent manner [Ref. 9]. This is the case of a directionally unstable ship 
which is well known in the literature [Ref. 3]. If the condition in (14) is not met, it means 
that there is at least one complex conjugate root with real parts and will result in an 
oscillatory response, indicating loss of stability. 
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To analyze the limiting case of loss of stability equation (14) must be solved in the 



following form: 



BCD-AD^ -EB^ =0 



(16) 



The result of this equation will be a curve of zg as a function of xg and will be our locus 
for loss of stability. We can express the coefficients of equation (16) in the form: 



Au Az, and As are of the form given in [Ref. 13, pg. 14]. 

B\ and £3 are of the form given in [Ref. 13, pg. 14]. With the addition of a pitch angle, w, 
B 2 takes the form: 



B, =-m(K,U\K -N,)-nAlu 

+ niYp ^ - Umx ^ )+ mK^ ~ ) 

-Um)+m(N^u\mx^ -Y,) 

-m{K,- 1 ^ lN,U)+m UYp {mx^ - N, ) 

-mU{Np+ Xm - Y , )+ mUK, {mx^ - N, ) 

+ mzc w(m -Y,\l ^)- mZc w{mXc - Y^ Xf^c ~^i) 






(17) 



B — B^Zg + ^2 ^3 



(18) 



C — Cl Z q ~^C 2 Zq + C3 



(19) 
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Cl and C 3 are of the form given in [Ref. 13, pg. 15]. With the addition of a pitch angle, 
w, C 2 takes the form: 



=mU{rm^ - N,lY^u)-mUYp{N^U)-mUK,{N,U) 

■^mU{N. +l^XY,U)-mU{N^u\m-Y,) 

+ W{m-Y,ll^-N,)-W{mx^-Y,Xmx^-N,) 

-m{X sB-x'cWXmXc -Y^)+m{UN^ -UmXcXXpU) 

-m{UY^ -UmXNpU)+m{K^UXUN^ -UmXc) 

-mZaw{m-Y,XUN , -Umx^)-mzMY,UXl ^ ~N,) 

-mzcw{mx(; -Yf X^v^ )+mZcw{mXc -N- XUY^ -Urn) 

D = DjZc+D2 (20) 

Di is of the form given in [Ref. 13, pg. 16]. With the addition of a pitch angle, w, Dz 
takes the form: 

D, =UK,(x,B-x^Wlm-Y,)-UKXN.u1y,u) 

^UK,{n ,u\y,U)-(K,-I ^lx,B-x„wiY,u) 

+ K,(x,B-x„WIuY, -U m)-(K.Ulx,B-x^Wlmx^ -Y,) 

+ (kMUN, -Umx^iY,u)-{K,uXuY, -VnifN fl) 

-(k^u\y,uIVN, -V mx^)+{K^U\UK, -UmlN.U) 

+ mZcw{Y,UXUN, -Umxc )-mz,w{UK, -UmXN.U) ' 

The equation for the coefficient, E, remains unchanged [Ref. 13, pg. 16]. Applying the 
stability criterion, equation (16), and utilizing them in the resulting fifth order 
polynomial, F, [Ref. 13, pg. 18] we are able to solve F using the MATLAB program in 
Appendix A. The curves show zo as a function of xq and we show results for varying 
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pitch angles, w, and varying forward velocities, U. On all of the following graphs the 
pitch angle is varied from 10 degrees to -10 degrees in increments of five degrees. The 
top curve is the 10 degree curve and the bottom curve is the -10 degree curve. 
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ZG in ft ZG in ft 




Figure 1: ZG vs. XG for U = 2 ft/sec 




Figure 2: ZG vs. XG for U = 3.5 ft/sec 
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ZG in i\ ZG in ft 




Figure 3: ZG vs. XG for U = 5 ft/sec 




Figure 4: ZG vs. XG for U = 6.5 ft/sec 
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Figure 5: ZG vs. XG for U = 8 ft/sec 
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From the results of Figures 1 through 5, the following conclusions can be drawn: 

1. In all cases, a sufficiently high metacentric height is required in order to ensure 
vehicle stability. Regions of parameters that fall below the critical curves 
correspond to dynamic instability. 

2. The critical metacentric height that is required for dynamic stability is an increasing 
function of both vehicle speed and trim angle. 

3. Static stability alone does not necessarily ensure dynamic stability of motion during 
the turn. 

4. The loss of stability experienced here is a dynamic loss of stability. At the critical 
metacentric height, one pair of complex conjugate eigenvalues possesses a zero real 
part. This is an oscillatory loss of stability, which can not be predicted by considering 
the uncoupled sway/yaw equations of motion alone. 
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IV. NONLINEAR ANALYSIS 



A. INTRODUCTION 

In the previous chapter we performed a linear analysis of the equations of motion, 
observing that with changes to specific parameters (xq, zg^ 'v) h is possible to pass 
through a stable region to a region of loss of stability. 

In previous work [Ref 13] it was shown that these bifurcations to periodic 
solutions were all supercritical. This means that limit cycles were produced after the loss 
of stability. By introducing a pitch angle, w, in the nonlinear analysis we will analyze 
whether the bifurcations to periodic solutions will remain supercritical and what changes 
occur to the limit cycles themselves. 

B. TfflRD ORDER EXPANSIONS 

Our linearized system was written in the form of equation (1 1), ignoring the 
nonlinear terms. Including the nonlinear terms changes the form of equation (1 1) to: . 

Ax = B'x + g(x) 



where: 



g(x)= 



'SiixJ 
giix) 
8 six) 



Keeping terms up to third order, the matrix g(x) can be written in the vector form: 

5 (^)+ g (^)+ c(x) (22) 

where contains the second order nonlinear terms and contains the third 

order nonlinear terms. The cross flow integrals and the second order nonlinear terms 
remain unchanged with the addition of a pitch angle, w. [Ref. 13, pg.22, 23] However 
the third order nonlinear terms take the form: 

8?\x) 

8?ix 
8?{x) 

gp) = -/p) --(w - B)l>^ cosd 

6 

8?=-I?-Ux^W-x,b)i>^ cose 

O 

g?=Uz,W-ZsBy cose 
o 

=0 

The Taylor series expansion yielding the second and third order linear terms of the cross 
flow integrals, and the inverse of the system matrix, {A')~' , remain unchanged. [Ref. 13, 



where: 



g(x)= 



20 



pg. 24, 25]. However the B’ matrix has changed as shown in the Linear Analysis section 
of this work. These changes to the 4 x 4 F matrix are: 





^\2 


E}L 




D 


D 


D 


D 
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r 22 


^23 


Ll 


D 


D 


D 


D 






f 

^ 33 


^34 


D 


D 


D 


D 


0 


0 


1 


0 



where; 



Fn=“,,(Y,U)*a,,(N,U)+a„(K.U) 

F\2 -mx^^J *a„(K, +mz„)j 

Fn =^n^,u)*a,,(N^Uya„(K^u) 

F ,4 = a,, (W —B)cosG +a^ 2 {x(;W — XgB)cos6 + aij{—Z(;W + ZbB)cos9 
F„ =a,,{Y.U)+a,,{N^U)+ajK,U) 

^22 — ^21 (^v ^22 ("^v )t/ ^23 v 

Fj, =‘‘2,(ypU]+a,,{N^u)+a^(K^u) 

F24 =«2i (W -F)cos 0 + 022(^0^ — X gB)cosQ + a 22 {— ZqW + ZbB)cosO 

F„ =a„(}'.[/)+«„(W.C/)+a„{A:.t/) 

F„ =a„{Y,-mp + a,2{N,-mx^'P + a„(K,+mz^yj 

F„ =ajY,u)+a,,(N,u)+a„(K^u) 

F34 = Oj, (W — F)cos 0 +a 32 (xglT — Xg 5 )cos 0 + a33(— Z(jiy + ZbB)cos 6 
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The remaining elements of the nonlinear and constant terms remain unchanged [Ref . 13, 
pg. 28, 29]. 

C. COORDINATE TRANSFORMATIONS 

To continue our analysis it is necessary to bring our transform our coordinate 
system from state space to a normal coordinate system. This transformation is performed 
in the manner given [Ref. 13, pg. 29, 30]. 



D. CENTER MANIFOLD EXPANSIONS 

The center manifold expressions are of the form given [Ref. 13, pg. 30-33]. 
HoNvever, the coefficients 2, 3; j=5, 6, 7 are: 



r ^11 { 2 2 \ 

^ 1,5 “ \^31 ^21 j 

+ +^yzW3,m2, +ycm2,m„) 



a 



13 



D 






31 






^1.6 =-^>'g(2W3i" 132 +2^2, -|-m22) 



^12 

D 



a 



2/x>"*3i'«32 +^yz ("*3l'«22 ) 

+ yci^2i^u+f^22^n) 

^;o('«3l"'22 +"J32"*21 yz^2^f^^2 

^ [+myc(7n„m32 +m^2m^^)+2{ycW - 



(23) 



(24) 



22 



ki =-^>’cte +"122) 



D 



+ +k’Z^^2^22+yc^22f^l2) 

-^t;o."^32"^22 +kzf^22 +mycm^^m^^ +{yc^ ~ yh^hk] 



k.5 = ^yG{fnu+ml,) 

+ — — (/ mji yz^2\^2\ "*"!Vc^21^1i) 



D 



( 



^23 

D 



^jc>"*3l'”2l -^kz^2\ +"«>’c'”ll'”31 



+ y^Wml^ -yeBml^ 



k.6 =-^>’c(2"l3l'«32 +2^21 +m22) 



^22 

~D 



‘23 



2/,^,m3,m32+/^,(, 

+ Jc("*2l'”l2+"^22"^ll) 



kV"^3i'«22+W32"*21, 



^:o("^3l"^22 +"*32'«21 )+2/j.,m2,m 



>2 



^22 



^ [+m>’c(wi,m32+m,2m2i)+2(j(-W-3;g5)m4,m42_ 
^2,7 =-^}'Gk32+"^22) 



+ ^ (^;r>.'«32 + kz^22^22 + yG"^22"^12 ) 

-•^tx>"^32"*22 +kz^l2 +my^m^^m^^+{y^W -ysB)ml^] 

k.5 =-^>'Gfe+'«2l) 

+ ^(^x>."^31 +^vz"^3l"^21 + yG"^2l"Jll ) 



( 25 ) 



(26) 



(27) 



(28) 
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(29) 



xy^ix^ix +Jyz^lx + + y aWml^ ~ y ^Bml) 



^3,6 ~ y G ^^' ix^il +2m2| +^22) 



a 



32 



D 



a 



33 

D 



2/,^m3,m32 +/3,,(m3,m22 +m32m2, ) 

+ }'g('”2i"^ 12 +"^22"^ll) 

^ X}X^3X^22 +^32^21 yz^2X^22 

+ "i>’c ('”ii'^32 +m,2m3, )+2(>cW - ygB)m^^m 



7 ^31 / 2 2 \ 

^3,7 “ G V^32 ^22 / 

xy^32 ^ yz^32^22 y G^22^\2^ 

-■ ^ t ; 0 .'” 32'«22 +^>:"^22 + "^>' g '« 12'«32 + (^ G ^^ “ > B ^ K 2 ] 



/ 



(30) 



(31) 



E. AVERAGING 

The procedure for averaging the equations up to the third order is the same one 
given in [Ref. 13, pg. 37-40]. The addition of a nonzero pitch angle yields new 
coefficients /,y, 7=1, 2, 3; y=l, 2, 3, 4. They are: 
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/ =-^ 
" D 



— —{EQinl^ +3E^mf^m2^ +3E2m^^ml^ +E2m\^^ 

6y 

+ -{ y / - B ) coseml , 

6 



a 



12 

D 



— +3£’2^u^2i +3£'3mjjm2i +£’4m2i) 

67 



’^ — { xoW - x ^ B ) cos9m 



oD 



(32) 



I _ ^11 

“ D 



67 



3 £o"^u"Ji 2 +3£] (^1^1^22 +2m,,m,2m2i ) 



+— (W — 5)cos03m4|m42 

6 



a 



12 

D 



+ 3£2 (^12^21 +2m2im22/«ii )+3-^3"^21^22 ; 

^3£imf|m,2 +3£'2(^u"J22 +2m,,m,2m2, )+3£#(m,2W2, +2m2im22»ii, ) 
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+ 3£4/n2,/W22 



+ — {xqW — Xg B)cos> 63ml^ 
6 



41^42 



+ ^ {zg^ - Zi,B)cosd3ml m ^2 



(33) 
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3£om,,ra,^2 +3£, (m,2m2, +2«2,,m,2Wi22)'*'^^2("^22"^ii +2m2, 0122 ^ 2 , 2 )'^ 



+ 3£3m2,m22 



+ — (W — 5)cos03m4,m 
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a 



12 



D 



C 



— {E^m]2 + 3 E 2 mf^m 2 ^ +3£3m,2m22 +£'4^22)+ — -XgB)cosBm\. 

6y 6 



^^i^GW-ZsB)coseml2 

bu 
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a 



21 

D 



■ — +3£,mf,m2, +3E2W,,m2, +E3m2, )+ — (W’-E)cosftn4 

67 6 

+3E2m,^,m2, +3E3m,,m2, +£ 4 / 722 , )+ — -JCg£)cosft72; 

67 6 



C 
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3Eo/ 72,^,/72,2 +3E, (/72,^,/T222 + 2/72, , 772, 2/2221 )+ 3Ej (/72,2/722, + 2/722, 77222^, , ) 
+ 3E3 7722 , 772 22 
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1 



— (w — E)cOS 037724 , 77242 

6 

3E,/72,^, 772,2 + 3 E 2 (/72,^, 77222 + 2/72, , 772,2/722, )+ 3 E 3 (/72,2/722, +2/72 
+ 3E4/722 , 772 22 



'2, I ^a/2 2, 772 22 772, , ^ 






— — XgE)cOS 03 / 724 , 

6 

jE)cos 03772; 



14,77/42 



14,77242 



(36) 



'^3Eo772,,772,2 + 3E, (t 72, 27722, + 2/n, , 772,277222 )+ 3E2 (77222772, , + 2/722,771 



12 , 77222 772,2 )Y 



+ 3E3 77/2 , 772 22 



(W^ — E)cos 0 377/4, 77/42 

/ 

3E,772,, 772,\ +3E2 (772,^2"^21 + 2/72,, 772,2/7222 )+ 3 E 3 (77222/72, , + 2/722,77222772,2 )^ 

_ O 



f 3 E 4 77/2 , 772 22 



IV —XgB)cos03m^^ml2 



+ — (zc^ “ Za 5 )cos 03 m 4 , 



D 



m 



42 
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a 



21 

D 



C 



— +3£,m,^,m2, +?>E 2 m^^ml 2 + — -B)cosdm, 

67 6 



a 



22 



Z) 



c. 
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(£,mf 2 + 3 £ 2 ^n^ 2 i + 'iE^m^ 2 ^]i +-^4^22) 



+ — (x^ W — Xg 5 )cos 6^4 
6 



+ -g^(z(jiy — Zg 5 )cos ^ 



( 38 ) 



/ — ^31 

Z) 



67 



(^oWif, + 3E^mf^m2^ +3£2^n^2i +-£'3W32i )+ — -5)cos6fn4, 

6 



a 



32 



z> 



— + 3 £' 2 mf,m 2 , + 3 £’ 3 m,,m 2 , +£4^2,) 

67 

+ — (x(j W-Xgfi)cos ftn4, 



+ — {zcW-ZsB)cosdm 



6D 



d -3 



( 39 ) 



/ — 31 

Z) 



Cr 



67 



3EQmf^m^2 +3£',(mf,m22 +2mi,m,2m2i )+3£2 (^ 12^21 + 2m2,m22W,, ^ 



+ 3 E 2 m\^m 22 



+—{W — B^cosd3m\m^2 



a 



32 

Z) 



c. 



67 



1 



3£,mf|7n,2 +3£2(wf,m22 +2m,,m,2m2, )+3£3(m,2m2i +2m2,m22"3n ) 



+ 3 E^m 2 im 22 



+ -{x(jW-XgB)cosd 3 m"^ 

6 



+ ^{zc^ -ZBB)coseml^m ^2 

bu 



( 40 ) 
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/ — ^31 

” ” D 



Cr 



6y 



+ 3E2 (^22^11 ■*■ +2m2i ^22 W[2 )+ 3£'3m2im 



2 

22 



a 



32 



D 



+ — (W — 5)cos03m4,m42 

6 

+3£2te"^2i +2mi,m,2m22) 

[+3E2{ml2m^^ +2m2^m22m^2)^3E^m2^ml2 

+ — {x^W — XgB)cos 63m^jm^2 

6 



■'■^(^c^-Z5^)cos03m4,m42 

ou 



(41) 



/ — ^31 

D 



C 



— (£'omf2 +3£'im,^jm2i +3£2^i2^22 ■*'■^ 3 ^ 22 ) 
67 

+ — (ly - 5 )cos 0m 42 

6 



+ ■ 



a 



32 



D 



C 



+ 3£2"^H^21 ■*■ '3E2^m^2^22 ■•■ ^4^22 ) 
67 

+ — - XgB)cos0ml2 
6 



+ -^(zrW - ZoB)cos6m 

6 ^ V G B ; 



(42) 



The remaining procedure for determining the equation for the radius of the resulting limit 
cycles is identical to the one described in [Ref. 13, pg. 43, 44], which is: 

R = a'eR + KR^ (43) 
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F. LIMIT CYCLE ANALYSIS 



At steady state, R = 0, equation (43) becomes: 

Q = R{a'e + KR^) 

Equation (66) has two solutions: 



R = 0 




(45) 

(46) 



Equation (45) is the trivial solution and gives no useful information. Equation (46) gives 
a constant amplitude limit cycle, R, in the cartesian coordinate system. This limit cycle 
will exist if the quotient inside the radical sign is positive: 



— a s 
K 



>0 



(47) 



This condition is necessary for /? to be a real number. Since a ’ is always negative, the 
existence of limit cycles depends on the parameter K. The introduction of a nonzero pitch 
angle does not change the dependence the limit cycle has on the parameter K. Stated, this 
dependence is: 

• If ^<0, periodic solutions exist and they are stable. 

• If AT>0, periodic solutions exist and they are unstable. 
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G. RESULTS AND DISCUSSION 



Results for the stability parameter, K, are presented in Figures 6 through 15 
Figures 6 through 10 provide results for A” at a constant forward speed, U, and varying 
pitch angles. The pitch angles were varied from positive 10 degrees to negative 10 
degrees in 5 degree increments. In Figures 6 through 10, the bottom curve represents 
solutions for positive 10 degrees and the top curve represents negative 10 degrees. It is 
clear that all values of K are negative indicating they are stable solutions. Notice that 
decreasing pitch angles the solutions for K become less negative, indicating that while 
they are stable, these solutions are less stable than those at the higher pitch angles. In 
Figures 1 1 through 15, solutions for the stability parameter, K, are again represented, but 
with the pitch angle held constant and the forward speed, U, varied from 2 ft/sec to 8 
ft/sec in 1.5 ft/sec increments. The bottom curve in Figures 1 1 through 15 represents U = 
2ft/sec. As forward speed increases, the curves representing the stability parameter K 
tend to become more negative, but in a more pronounced fashion than those where U was 
held constant. 
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K*GAMMA K'GAMMA 



U=2 ft/sec 




Figure 6 : XG vs. K*GAMMA for U=2 ft/sec 



U=3.5 ft/sec 




Figure 7: XG vs. K*GAMMA for U=3.5 ft/sec 
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K'GAMMA K'GAMMA 



X 1 0”^ 




Figure 8: XG vs. K*GAMMA for U=5 ft/sec 



U=6. 5 ft/sec 




Figure 9: XG vs. K*GAMMA for U=6.5 ft/sec 
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2 



U=8 ft/sec 



xIO'^ 




Figure 10: XG vs. K*GAMMA for U=8 ft/sec 
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K-GAMMA K'GAMMA 



THETA=10deg 




Figure 11: XG vs. K*GAMMA for THETA =10 deg 



THETA=5 deg 




Figure 12: XG vs. K*GAMMA for THETA=5 deg 



34 



K'GAMMA K'GAMMA 



THETA=0 deg 



X 10'^ 




Figure 13: XG vs. K*GAMMA for THETA=0 deg 



THETA=-5 deg 




Figure 14: XG vs. K*GAMMA for THETA=-5 deg 
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THETA=-10deg 




Figure 15: XG vs. K*GAMMA for THETA=-10 deg 
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V. CONCLUSIONS AND RECOMMENDATIONS 



This thesis presented a continuing study of the formation of limit cycles due to the 
coupling of the sway/yaw/roll equations of motion. We have shown that loss of stability 
occurs in the form of stable limit cycles and that the addition of a nonzero pitch angle does 
not significantly affect the formation of these limit cycles. Through a linear analysis of the 
sway/yaw/roll equations of motion we demonstrated that the addition of a nonzero pitch 
angle affects the domain of stability of straight line motion, especially at higher speeds. 
This was validated by the nonlinear analysis as well. As a recommendation for further 
study in this area we suggest that the analysis be expanded to include coupling into the 
heave/pitch directions of motion as well as the effects automatic path control. 
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APPENDIX 



The following is a list of the MATLAB and FORTRAN programs used in this thesis. 
Complete printouts of the programs accompany this list. 

• THESIS1.M 

A MATLAB program for performing the linear analysis section of this thesis. 

• HOPF_lNEW.FOR 

A FORTRAN program for performing the nonlinear analysis section of this thesis. 
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% THESIS l.M 

% 

% LOSS OF STABILITY 

clear 

a=l; 

W= 12000; 

DCX=1760; IYY=9450; 

IZZ=10700;IXZ=0; IXY=0; 

IYZ=0; L=17.425; RHO=1.94; 

G=32.2; U=8.0; M=W/G; B=W; 

THETA=0; 

THETA=THETA*pi/l 80; 
OMEGA=U*tan(THETA); 

ND1=0.5*RHO*L''2; 

% DEFINE HYDRODYNAMIC COEFHCIENTS 

YPDOT= 1 .270e-04*ND 1 *L^2; 
YVDOT=-5.550e-02*ND 1 *L; 

YRDOT= 1 ,240e-03*ND 1 *L^2; 
YP=3.055e-03*NDl*L; 

YPOMEGA=0; 

YP=YP*U+YPOMEGA*OMEGA+M*OMEGA; 

YV=-9.310e-02*NDl; 

YVOMEGA=0; 

YV=YV+YV0MEGA*0MEGA; 

YR=-5.940e-02*NDl*L; 

YROMEGA=0; 

YR=YR+YROMEGA*OMEGA; 

NPDOT=-3.370e-05*NDl*L^3; 

N VDOT= 1 .240e-03 *ND 1 *L^2; 
NRDOT=-3.400e-03*ND 1 *L'^3; 
NP=-8.405e-04*ND 1 *L^2; 

NPOMEGA=0; 

NP=NP+N1^0MEG A* OMEGA; 
NV=-1.484e-02*NDl*L; 

NVOMEGA=0; 

NV=NV+NVOMEGA*OMEGA; 

NR=- 1 ,640e-02*ND 1 *L^^2; 

NROMEGA=0; 

NR=NR+NROMEGA*OMEGA; 
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KPDOT=- 1 .0 1 e-03*ND 1 
KVDOT= 1 ,27e-04*ND 1 *L^2; 

KRDOT=-3.37e-05*NDl *L^3; 

KP=- 1 . 1 0e-02*ND 1 *L'"2; 

KPOMEGA=0; 

KP=KP+KPOMEGA*OMEGA; 

KV=3.055e-03*NDl*L; 

KVOMEGA=0; 

KV=KV+KVOMEGA*OMEGA; 

KR=-8.41e-04*NDl*L^2; 

KROMEGA=0; 

KR=KR+KROMEGA*OMEGA; 

flag=0; 

for XG=0:0.01:2, 

flag=flag+l; 

xg(flag)=XG; 

a=IXX-KPDOT; b=KP*U; e=KV*U; 
f=KRDOT; i=YP*U; j=M-YVDOT; k=YV*U; 
l=XG*M-YRDOT; m=U*(YR-M); o=NPDOT; 
p=NP*U; q=-XG*W; r=XG*M-NVDOT; 
w=U*(NR-XG*M); x=NV*U; u=IZZ-NRDOT; 

al=-u*M'^2; 

a2=-u*M*YPDOT-u*M*KVDOT+M*o*l+f*r*M; 

a3=a*j*u-a*l*r-u*KVDOT*YPDOT+KVDOT*o*l+f*r*YPDOT-f*o*j; 

bl=(w*(M^2))+(r*U*(M^2)); 

b2=-M*e*u-M*u*i+w*M*YPDOT+w*KVDOT*M-o*m*M+p*l*M- 

f*x*M+r*M*U*YPDOT-o*j*M*U+r*KR*U*M; 

b3=-e*u*YPDOT+e*o*l-u*i*KVDOT+w*KVDOT*YPDOT- 

KVDOT*o*m+KVDOT*p*l-a*j*w-a*k*u+a*l*x+a*r*m-b*j*u+b*l*r+f*r*i- 

f*x*YPDOT+o*k*f-f*p*j+r*KR*U*YPDOT; 

b2=b2+M*OMEGA*j*u-M*OMEGA*l*r; 

cl=-x*(M'^2)*U; 

c2=r*i*M*U-x*M*U*YPDOT-x*KR*U*M+o*k*M*U-p*j*M*U+j*u*W-l*r*W- 

q*l*M+w*i*M-m*p*M+e*w*M; 

c3=a*k*w-a*m*x+b*j*w+b*k*u-b*l*x-b*r*m+r*i*KR*U- 

x*KR*U*YPDOT+o*k*KR*U-p*j*KR*U-q*l*KVDOT+w*i*KVDOT-m*p*KVDOT- 

e*u*i+e*w*YPDOT-e*o*m+e*p*]+f*q*j-f*x*i+f*p*k; 



41 



c2=c2-M*OMEGA*j*W-M*OMEGA*k*u-M*OMEGA*l*x+M*OMEGA*r*m; 



d2=q*j*M*U-x*i*M*U+p*k*M*U+q*m*M-j*w*W-k*u*W+l*x*W+r*m*W; 

d3=q*j*KR*U-x*i*KR*U+p*k*KR*U-f*q*k+q*m*KVD0T-e*q*l+e*w*i-e*m*p- 

b*k*w+b*m*x; 

d2=d2+M*OMEGA*k*w-M*OMEGA*U*(KR-M)*x; 

e2=k*w*W-m*x*W-q*k*M*U; 

e3=e*q*m-q*k*KR*U; 

f5=(-e2*(b 1 '^2))+b 1 *c 1 *d2; 

f4=((b2*c 1 +b 1 *c2)*d2)+b Pci *d3-a P(d2^2)-e3*(b 1 '^2)-2*e2*b 1 *b2; 
f3=-a2*(d2''2)-2*d3 =^d2*al -2*e3*b 1 *b2- 

e2*((b2^2)+2*b 1 *b3)+d2*(b3 *c l+b2*c2+b 1 *c3)+d3 *(b2*c 1 +b 1 *c2); 

f2=-e3*((b2^2)+2*bPb3)-2*e2*b2*b3+d2*(b3*c2+b2*c3)+d3*(b3=^cl+b2*c2+bPc3)- 

a3*(d2^2)-2*d3*d2*a2-aP(d3''2); 

fl=b3*c3*d2+d3*(b3*c2+b2*c3)-2*e3=*=b2*b3-e2*(b3''2)-2*d3*d2*a3-a2*(d3^2); 

fO=b3*c3*d3-a3*(d3'^2)-e3*(b3^2); 



coef=[f5 f4 f3 f2 fl fO]; 
ZG=roots(coef); 



tot(flag)=ZG(5,l); 

end 



42 



C PROGRAM HOPF_lNEW.FOR 

C EVALUATION OF HOPF BIFURCATION FORMULAS USING THE SUBOFF 
C SUBMARINE MODEL 
C 

IMPLICIT DOUBLE PRECISION (A-H,0-Z) 

REALMS L,IYY,M,YPDOT,YVDOT,NDl,YRDOT 
REALMS YP,YV,YR,NPDOT,NVDOT,NRDOT,NP,NV,NR 
REAL*8 GAMA,U,KVDOT,KRDOT,KPDOT,D 
REALMS E0,E1 ,E2,E3,E4,XG,ZG,KR,KP,KV,XG1 ,ZG1 
C 

REAL*8 Ml 1,M12,M13,M14,M21,M22,M23 
REAL*8 M24,M3 1,M32,M33,M34,M41,M42,M43,M44 
REAL*8 N1 1,N12,N13,N14,N21,N22,N23,N24 
REAL*8 N3 1 ,N32,N33,N34,N41 ,N42,N43,N44 
REAL*8 LI 1,L12,L13,L14,L21,L22,L23,L24,L31 
REAL*8L32,L33,L34,L15,L16,L17,L25,L26,L27,L35 
REAL*8L36,L37,L1A,L2A,L3A,L4A,L5A,L6A,L7A,L8A,L9A 
REAL*8 L10A,L1 1A,L12A,L1,L2,L3,L4,L5,L6,L7 
REAL*8 L8,L9,R1 1,R12,R13,R14,R21,R22,R23,R24 
REAL*8 PI 1,P12,P13,P21,P22,P23 
DOUBLE PRECISION THETA 
C 

DIMENSIONF(4,4),T(4,4),TINV(4,4),FVl(4),IVl(4),YYY(4,4) 

DIMENSION WR(4),WI(4),TSAVE(4,4),TLUD(4,4),IVLUD(4),SVLUD(4) 
DIMENSION AS AVE(4,4), A2(4,4),XL( 1 8),HT( 1 8),ZGG(202),FF(4,4) 
DIMENSION VEC0( 1 8),VEC 1 ( 1 8), VEC2( 1 8), VEC3( 1 8), VEC4( 1 8),XGG(202) 

C 

INTEGER U,K 
C 

OPEN (20,FILE='HOPF25.RES') 

. OPEN (2 1 ,FILE=’DAT25.DAT',STATUS='OLD') 

C OPEN (23,FILE='HOPFl.RES',STATUS='OLD') 

C 

WEIGHT= 12000.0 
DCX =1760.0 
lYY =9450.0 
IZZ =10700.0 
EXZ =0.0 
EXY =0.0 
lYZ =0.0 
L =17.425 
RHO =1.94 
CD =0.5 

C CD =0.5*CD*RHO ??? 

G =32.2 
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XB =0.0 
ZB =0.0 
YG =0.0 
YB =0.0 
YDELTAR=0.0 
DELTAR=0.0 
NDELTAR=0.0 
NPROP=0.0 
M =WEIGHT/G 
B = WEIGHT 
W=WEIGHT 

WRITE (*,*) ’ ENTER U' 

READ (*,*) U 

WRITE (*,*) ' ENTER THETA (DEGREES)' 
READ (*,*) THETA 
C U =5.0 

NDl = 0.5*RHO*L**2 
C THETA=5 

THETA=THETA*PF1 80 
OMEGA=U*DTAN(THETA) 

C 

C DEFINE HYDRODYNAMIC COEPHCIENTS 
YPDOT= 1 .270E-04*ND 1 
YVDOT=-5.550E-02*ND1 *L 
YRDOT= 1 .240E-03*NDI *L**2 
YP=3.055E-03*ND1*L 
YPOMEGA=O.DO 

YP=YP+YPOMEGA*OMEGA+M*OMEGA 

YV=-9.310E-02*ND1 

YVOMEGA=O.DO 

yv=yv+yvomega*omega 

YR=-5.940E-02*ND1*L 

YROMEGA=O.DO 

YR=YR+YROMEGA*OMEGA 

C 

NPDOT=-3.370E-05*ND 1 *3 

N VDOT= 1 .240E-03 *ND 1 *2 

NRDOT=-3.400E-03*ND 1 *L* *3 
NP=-8.405E-04*ND I *L**2 
NPOMEGA=O.DO 
NP=NP+NPOMEGA*OMEGA 
NV=-1.484E-02*ND1*L 
NVOMEGA=O.DO 

nv=nv+nvomega*omega 

NR=- 1 .640E-02*ND 1 *2 
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NROMEGA=O.DO 

NR=NR+NROMEGA*OMEGA 



C 



C 

C 

C 



C 



KPDOT=-1.01E-03*NDPL**3 

KVDOT= 1 ,27E-04*ND 1 

KRDOT=-3.37E-05*ND1*L**3 

KP=- 1 . 10E-02*ND 1 *L**2 

KPOMEGA=O.DO 

KP=KP+KPOMEGA*OMEGA 

KV=3.055E-03*ND1*L 

KVOMEGA=O.DO 

KV=KV+KVOMEGA*OMEGA 

KR=-8.4 1E-04*ND 1 *L**2 

KROMEGA=O.DO 

KR=KR+KROMEGA*OMEGA 

DEFINE THE LENGTH X AND HEIGHT H TERMS FOR THE INTEGRATION 

ALL IN FEET. 

XL( 1)=-105.9/12.0 
XL( 2)=- 104.3/ 12.0 
XL( 3)=-99.3/12.0 
XL( 4)=-94.3/12.0 
XL( 5)=-87.3/12.0 
XL( 6)=-76.8/12.0 
XL( 7)=-66.3/12.0 
XL( 8)=-55.8/12.0 
XL( 9)=72.7/12.0 
XL(10)=79.2/12.0 
XL(11)=83.2/12.0 
XL(12)=87.2/12.0 
XL(13)=91.2/12.0 
XL(14)=95.2/12.0 
XL(15)=99.2/12.0 
XL(16)=101.2/12.0 
XL(17)=102.1/12.0 
XL(18)=103.2/12.0 

HT(1)= 0.000 
HT( 2)= 2.28/12.0 
HT( 3)= 8.24/12.0 
HT( 4)= 13.96/12.0 
HT( 5)= 19.76/12.0 
HT( 6)= 25.1/12.0 
HT( 7)= 29.36/12.0 
HT( 8)= 31.85/12.0 
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HT( 9)= 31.85/12.0 
HT(10)= 30.00/12.0 
HT(11)= 27.84/12.0 
HT(12)= 25.12/12.0 
HT( 13)= 21.44/12.0 
HT(14)= 17.12/12.0 
HT(15)= 12.0/12.0 
HT(16)= 9.12/12 
HT(17)= 6.72/12 
HT(18)= 0.00 
C 

DO 104 K= 1,18 
VEC0(K)=HT(K) 

VEC 1 (K)=XL(K)*HT(K) 
VEC2(K)=XL(K)*XL(K)*HT(K) 
VEC3(K)=XL(K)*XL(K)*XL(K)*HT(K) 
VEC4(K)=XL(K)*XL(K)*XL(K)*XL(K)*HT(K) 
104 CONTINUE 
CALL TRAP(1 8,VEC0,XL,E0) 

CALL TRAP( 1 8, VEC 1 ,XL,E 1 ) 

CALL TRAP( 1 8, VEC2,XL,E2) 

CALL TRAP(18,VEC3,XL,E3) 

CALL TRAP( 1 8, VEC4,XL,E4) 

C 

GAMA=0.001 



C READ THE CRIHC AL VALUES FOR XG AND ZG FROM FILE DATA 1 .DAT 
C 

XGG(1)=0.0 
ZGG(1)=0.01 6358083 
DO 1 1=2,202 
READ (21,*)XG,ZG 
C WRITE(*,*)XG,ZG 
XGG(I)=XG 
ZGG(I)=ZG 



C DETERMINE [F] COEFFICIENTS 
C 

D=((M-YVDOT)*(IZZ-NRDOT)*(IXX-KPDOT)) 

&-((M-YVD0T)*(DCZ-KRD0T)*(-IXZ-NPD0T» 

&-((M*XG-NVD0T)*(M*XG-YRD0T)*(IXX-KPD0T)) 

&+((M*XG-NVDOT)*(IXZ-KRDOT)*(-M*ZG-YPDOT)) 

&+((-M*ZG-KVDOT)*(M*XG-YRDOT)*(-IXZ-NPDOT)) 
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&-((-M*ZG-KVD0T)*(IZZ-NRD0T)*(-M*ZG-YPD0T)) 

c 

A 1 1 =((IZZ-NRDOT)*(IXX-KPDOT))-((IXZ-KRDOT)*(-IXZ-NPDOT)) 

A12=((-M*XG+YRD0T)*(IXX-KPD0T))+((IXZ-KRD0T)*(-M*ZG-YPD0T)) 

A13=((M*XG-YRD0T)*(-]XZ-NPD0T))-((IZZ-NRD0T)*(-M*ZG-YPD0T)) 

A2 1 =((-M*XG+NVDOT)*(IXX-KPDOT))+((-M*ZG-KVDOT)*(-IXZ-NPDOT)) 

A22=((M-YVDOT)*(IXX-KPDOT))-((-M*ZG-KVDOT)*(-M*ZG-YPDOT)) 

A23=((-M+YVDOT)*(-IXZ-NPDOT))+((M*XG-NVDOT)*(-M*ZG-YPDOT)) 

A3 1 =((M*XG-NVDOT)*(IXZ-KRDOT))-((-M*ZG-KVDOT)*(IZZ-NRDOT)) 
A32=((-M+YVDOT)*(IXZ-KRDOT))+((-M*ZG-KVDOT)*(M*XG-YRDOT)) 
A33=((M-YVDOT)*(IZZ-NRDOT))-((M*XG-NVDOT)*(M*XG-YRDOT)) 
C============================================================== 

C EVALUATE TRANSFORMATION MATRIX OF EIGENVECTORS 

C 

F( 1 , 1 )=(A 1 1 * YV*U+A 1 2*NV*U+A 1 3 *KV*U)/D 

F(1,2)=(A11*(YR*U-M*U)+A12*(-M*XG*U+NR*U)+A13*(M*ZG*U+KR*U))/D 
F( 1 ,3)=( A 1 1 * YP*U+A 1 2*NP*U+A 1 3* (KP*U-M*ZG*OMEG A))/D 
F(l,4)=(All*(W-B)*DCOS(THETA)+A12*XG*(W-XB)*B*DCOS(THETA)+ 
&A13*(-ZG*W+ZB)*B*DCOS(THETA))/D 
F(2, 1 )=( A2 1 * YV*U+A22*NV*U+A23*KV*U)/D 

F(2,2)=(A21*(YR*U-M*U)+A22*(-M*XG*U+NR*U)+A23*(M*ZG*U+KR*U))/D 

F(2,3)=(A2i*YP*U+A22*NP*U+A23*(KP*U-M*ZG*OMEGA))/D 

F(2,4)=(A21*(W-B)*DCOS(THETA)+A22*(XG*W-XB*B)*DCOS(THETA)+ 

&A23*(-ZG*W+ZB*B)*DCOS(THETA))/D 

F(3, 1 )=(A3 1 * YV*U+A32*NV*U+A33*KV*U)/D 

F(3,2)=(A3 1 *(YR*U-M*U)+A32*(-M*XG*U+NR*U)+A33*(M*ZG*U+KR*U))/D 
F(3,3)=(A3 1 * YP*U+A32*NP*U+A33*(KP*U-M*ZG*OMEGA))/D 
F(3,4)=(A3 1 *(W-B)*DCOS(THETA)+A32*(XG*W-XB=*=B)*DCOS(THETA)+ 
&A33*(-ZG*W+ZB*B)*DCOS(THETA))/D 
F(4,l)=0.0 
F(4,2)=0.0 
F(4,3)=1.0 
F(4,4)=0.0 
C 

D0 11K=1,4 
DO 12 J=l,4 
ASAVE(K,J)=F(K,J) 

12 CONTINUE 
11 CONTINUE 

CALL RG(4,4,F, WR,WI, 1 , Y YY,IV 1 ,FV 1 ,IERR) 

WR]TE(23, 1 007)WR( 1 ),WR(2),WR(3),WR(4) 

CALL DSOMEG(IEV,WR,WI, OMEGA, CHECK) 

C WRITE (*,*) CHECK 

C WRITE (*,*) WR(2) 
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C WRITE (*,*) WI(4) 

OMEGAO=OMEGA 
D0 5 J=l,4 
T(J,1)= YYY(J,ffiV) 

T(J,2)=-YYY(J,IEV+1) 

5 CONTINUE 

IF (IEV.EQ.1.0) GO TO 13 
IF (IEV.EQ.2.0) GO TO 18 
IF (IEV.EQ.3.0) GO TO 14 
C STOP 3004 
14 D0 6J=1,4 
T(J,3)=YYY(J,1) 

T(J,4)=YYY(J,2) 

6 CONTINUE 
GO TO 17 

18 DO 19 J=l,4 
T(J,3)=YYY(J,1) 

T(J,4)=YYY(J,4) 

19 CONTINUE 
GO TO 17 

13 D0 16J=1,4 
T(J,3)=YYY(J,3) 

T(J,4)=YYY(J,4) 

16 CONTINUE 

17 CONTINUE 

NORMALIZATION OF THE CRITICAL EIGENVECTOR 
CALL NORMAL(T) 

INVERT TRANSFORMATION MATRIX 

D0 2K=1,4 
DO 3 J=l,4 
TINV(K,J)=0.0 
TSAVE(K,J)=T(K,J) 

3 CONTINUE 
2 CONTINUE 

CALLpLUD(4,4,TSAVE,4,TLUD,IVLUD) 

C DO 4 J= 1,4 

C IF (IVLUD(J).EQ.O) STOP 3003 
C 4 CONTINUE 

CALL DILU(4,4,TLUD,IVLUD,SVLUD) 

D0 8K=1,4 
D0 9 J=l,4 
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TINV(K,J)=TLUD(K,J) 

9 CONTINUE 
8 CONTINUE 

CHECK Inv(T)*A*T 

CALL MULT(TINV,AS AVE,T,A2) 

P1=A2(3,3) 

P2=A2(4,4) 

PEIG1=P1 

PEIG2=P2 

DEFINITION OF Nij 

N11=TINV(1,1) 

N12=TINV(1,2) 

N13=TINV(1,3) 

N14=TINV(1,4) 

N21=TINV(2,1) 

N22=TINV(2,2) 

N23=TINV(2,3) 

N24=TINV(2,4) 

N31=TINV(3,1) 

N32=TINV(3,2) 

N33=TINV(3,3) 

N34=TINV(3,4) 

N41=TINV(4,1) 

N42=TINV(4,2) 

N43=TINV(4,3) 

N44=TINV(4,4) 

DEFINITION OF Mij 

M11=T(1,1) 

M12=T(1,2) 

M13=T(1,3) 

M14=T(1,4) 

M21=T(2,1) 

M22=T(2,2) 

M23=T(2,3) 

M24=T(2,4) 

M31=T(3,1) 

M32=T(3,2) 

M33=T(3,3) 
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M34=T(3,4) 

M41=T(4,1) 

M42=T(4,2) 

M43=T(4,3) 

M44=T(4,4) 

C 

C DEFINITION OF Lij 
C 

L1=YG*((M3 P*2)+(M2 1 * *2)) 

L2=DCY*(M3 1 **2)+IYZ*M3 1 *M2 1 +YG*M2 1 *M 1 1 
L3=DCY*M31*M21+IYZ*(M21**2)+M*YG*M11*M31+YG*W*(M41**2)- 
YB*B*(M41** 

& 2 ) 

L4=YG*(2.0*M3 1 *M32+2.0*M2 1 *M22) 

L5=2.0*IX Y*M3 1 *M32+IYZ*(M3 1 *M22+M32*M2 1)+YG*(M2 1 *M 1 2+M22*M 1 1 ) 

L6=KY*(M3 1 *M22+M32*M21 )+2.0*IYZ*M21 *M22+M* YG*(M 1 1 *M32+M 12*M3 1 
)+ 2 . 

&0*(YG*W-YB*B)*M41*M42 

L7=YG*((M32**2)+(M22**2)) 

L8=IXY*(M32**2)+IYZ*M32*M22+YG*M22*M12 

L9=IXY*M32*M22+IYZ*(M22**2)+M*YG*M12*M32+(YG*W- 

YB*B)*(M42**2) 

C 

C 

L15=(A11/D)*L1+(A12/D)*L2-(A13/D)*L3 
L 1 6=( A 1 l/D)*L4+( A1 2/D)*L5-(A 1 3/D)*L6 
LI 7=( A 1 l/D)*L7+( A 1 2/D)*L8-(A 1 3/D)*L9 
L25=(A21/D)*L1+(A22/D)*L2-(A23/D)*L3 
L26=(A21/D)*L4+(A22/D)*L5-(A23/D)*L6 
L27=(A21/D)*L7+(A22/D)*L8-(A23/D)*L9 
L35=(A31/D)*L1+(A32/D)*L2-(A33/D)*L3 
L36=(A31/D)*L4+(A32/D)*L5-(A33/D)*L6 
L37=(A31/D)*L7+(A32/D)*L8-(A33/D)*L9 
C 

C=CD/(6.0*GAMA) 

LI A=C*(E0*(M 1 1 **3)+3.0*E 1 *(M 1 1 **2)*M2 1+3.0*E2*M 1 1 *(M2 1 **2)+E3*(M2 1 
&**3))-((1.0/6.0)*(W-B)*DCOS(THETA))*(M41**3) 

L2A=C*(E 1 *(M 1 1 **3)+3.0*E2*(M 1 1 **2)*M2 1+3.0*E3*M 1 1 *(M2 1 **2)+E4*(M2 1 
&**3))-((1.0/6.0)*(XG*W-XB*B)*DCOS(THETA))*(M41**3) 
L3A=((1.0/6.0)*(ZG*W-ZB*B)*DCOS(THETA))*(M41**3) 
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L4A=C*(3.0*E0*(M 1 1 **2)*M 1 2+3.0*El *((M 1 1 **2)*M22+2*M 1 1 *M 1 2*M2 1 )+3.0 
&*E2*(M12*(M21**2)+2.0*M21*M22*Mll)+3.0*E3*(M21**2)*M22)-((1.0/6.0) 
&*(W-B)*DCOS(THETA))*3.0*(M41**2)*M42 

L5 A=C*(3.0*E 1 *(M 1 1 **2)*M 1 2+3.0*E2*((M 1 1 **2)*M22+2*M 1 1 *M 1 2*M2 1 )+3.0 
&*E3*(M 12*(M2 1 **2)+2*M21 *M22*M1 1 )+3.0*E4*(M21 **2)*M22)-((1 ,0/6.0)*(X 
&G*W-XB*B)*DCOS(THETA))*3.0*(M41**2)*M42 
L6A=((1.0/6.0)*(ZG*W-ZB*B)*DCOS(THETA))*3.0*(M41**2)*M42 

L7A=C*(3*E0*M11*(M12**2)+3.0*E1*((M12**2)*M21+2*M11*M12*M22)+3.0*E 

&2*((M22**2)*M11+2.0*M21*M22*M12)+3.0*E3*M21*(M22**2)H(1.0/6.0)*(W 

&-B)*DCOS(THETA))*3.0*M41*(M42**2) 

L8A=C*(3.0*El*Mll*(Mi2**2)+3.0*E2*((M12**2)*M21+2.0*Mll*M12*M22)+3 
&.0*E3*((M22**2)*M11+2*M21*M22*M12)+3.0*E4*M21*(M22**2)) 
L9A=((1.0/6.0)*(ZG*W-ZB*B)*DCOS(THETA))*3.0*M41*(M42**2) 
L10A=-C*(E0*(M12**3)+3.0*E1*(M11**2)*M21+3.0*E2*M12*(M22**2)+E3*(M 
&22**3)H(1 ,0/6.0)*(W-B)*DCOS(THETA))*(M42**3) 

L11A=-C*(E1*(M12**3)+3.0*E2*(M11**2)*M21+3.0*E3*M12*(M22**2)+E4*(M 

&22**3)H(1.0/6.0)*(XG*W-XB*B)*DCOS(THETA))*(M42**3) 

L12A=((1.0/6.0)*(ZG*W-ZB*B)*DCOS(THETA))*(M42**3) 

C 

LI 1=(-A1 1/D)*L1 A+(A12/D)*L2A+(A13/D)*L3A 
LI 2=(- A 1 1/D)*L4A+(A 1 2/D)*L5 A+(A 1 3/D)*L6A 
L13=(-A11/D)*L7A+(A12/D)*L8A+(A13/D)*L9A 
L14=(-A1 1/D)*L10A+(A12/D)*L1 1A+(A13/D)*L12A 
C 

L2 1=(-A2 1/D)*L1 A+(A22/D)*L2A+(A23/D)*L3 A 
L22=(- A2 1/D)*L4A+(A22/D)*L5 A+(A23/D)*L6 A 
L23=(-A21/D)*L7A+(A22/D)*L8A+(A23/D)*L9A 
L24=(-A21/D)*L10A+(A22/D)*L11A+(A23/D)*L12A 
C 

L3 1=(-A3 1/D)*L1 A+(A32/D)*L2A+(A33/D)*L3A 
L32=(-A31/D)*L4A+(A32/D)*L5A+(A33/D)*L6A 
L33=(-A31/D)*L7A+(A32/D)*L8A+(A33/D)*L9A 
L34=(- A3 1 /D)*L1 0A+(A32/D)*L1 1 A+(A33/D)*L 1 2 A 
C 

R1 1=(N1 1*L1 1)+(N12*L21)+(N13*L31) 

R 1 2=(N 1 1 *L1 2)+(N 1 2*L22)+(N 1 3*L32) 

R 1 3=(N 1 1 *L 1 3)+(N 1 2*L23)+(N 1 3*L33) 

R 1 4=(N 1 1 *L1 4)+(N 1 2*L24)+(N 1 3*L34) 

R2 1 =(N2 1 *L1 1 )+(N22*L2 1 )+(N23*L3 1 ) 

R22=(N2 1 *L1 2)+(N22*L22)+(N23*L32) 

R23=(N2 1 *L1 3)+(N22*L23)+(N23*L33) 
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R24=(N21*L14)+(N22*L24)+(N23*L34) 

C 

P11=(N11*L15)+(N12*L25)+(N13*L35) 

P 1 2=(N 1 1 *L 1 6)+(N 1 2*L26)+(N 1 3 *L36) 

P 1 3=(N 1 1 *L 1 7)+(N 1 2*L27)+(N 1 3 *L37) 

P2 1=(N2 1 *L 1 5)+(N22*L25)+(N23 *L35) 

P22=(N2 1 *L1 6)+(N22*L26)+(N23*L36) 

P23=(N2 1 *L 1 7)+(N22*L27)+(N23 *L37) 

EVALUATE DALPHA AND DOMEGA 

ZGR =ZGG(I) 

ZGL=ZGG(I-1) 

ZGl =ZGR 
XGl =XGG(I) 

CALL FM ATRIX(ZG 1 ,XG 1 ,FF,U, THETA) 

CALL RG(4,4,FF,WR,WI,0, YYY,IV 1 ,FV 1 ,IERR) 

CALL DSTABL(DEOS,WR,WI, FREQ) 

ALPHR=DEOS 
OMEGR=FREQ 

ZGl =ZGL 
XGl =XGG(I-1) 

CALL FM ATRK(ZG 1 ,XG 1 ,FF,U, THETA) 

CALL RG(4,4,FF,WR,WI,0, YYY,IV 1 ,FV 1 ,ffiRR) 

CALL DSTABL(DEOS,WR,WI,FREQ) 

ALPHL=DEOS 
OMEGL=FREQ 

DALPHA=(ALPHR-ALPHL)/(ZGR-ZGL) 
DOMEGA=(OMEGR-OMEGL)/(ZGR-ZGL) 

EVALUATION OF HOPF BIFURCATION COEITICIENTS 

COEF 1 =( 1 .0/8 .0) * (3 ,0*R 1 1 +R 1 3+R22+3 ,0*R24) 
COEF2=(1.0/8.0)*(3.0*R11+R23-R12-3.0*R14) 

C AMU2 =-COEFl/(8.0*DALPHA) ???????? 

C BETA2=0.25*COEF1 ??????? 

C TAU2 =-(COEF2-DOMEGA*COEFl/DALPHA)/(8.0*OMEGAO) 
C PER =2.0*3. 1415927/OMEGAO 
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C PER =PER*U/L 

WRITE (20, 200DXG, ZG, COEFl ,D ALPHA, OMEGAO,PEIG1 ,PEIG2 
C WRITE (20,200 DCOEFl 

1 CONTINUE 
C 

STOP 

2001 FORMAT (7E14.5) 

4001 FORMAT (3F15.5) 

1007 FORMAT (4E14.5) 

END 
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